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Abstract 

Many physical and biological systems exhibit intrinsic cyclic dynamics that are altered by 
random external perturbations. We examine continuous-time autonomous dynamical systems 
exhibiting a stable limit cycle, perturbed by additive Gaussian white noise. We derive a formal 
approximation for the dynamics of sample paths that stay close to the limit cycle, in terms of 
a phase coordinate and a deviation perpendicular to the limit cycle. To leading order in the 
deviation, the phase advances at the deterministic speed superimposed by a Brownian-motion- 
like drift. The deviation itself takes the form of an (n — l)-dimensional Ornstein-Uhlenbeck 
process. We apply these results to the case of limit cycles emerging through a supercritical Hopf 
bifurcation, which is widespread in ecological and epidemiological models. We derive approx¬ 
imation formulas for the system’s stationary autocovariance and power spectral density. The 
latter two reflect the effects of perturbations on the temporal coherence and spectral bandwidth 
of perturbed limit cycles. We verify our results using numerical simulations and exemplify their 
application to the El Nino Southern Oscillation. 
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1 Introduction 

Many deterministic mathematical models exhibit stable limit cycles that are used to describe a 
wide spectrum of natural phenomena, ranging from population cycles [1-3], chemical oscillations 
[4], geophysical cycles [5] to periodic epidemic outbreaks [6, 7] and genetic oscillators [8]. Systems 
exhibiting intrinsically emerging oscillations are typically subject to external perturbations (e.g. 
the random change of an environmental parameter) or internal stochastic.ity (e.g. due to finite 
population sizes). Models are often extended to include such random factors by adding a stochastic 
term to the deterministic equations, typically in the form of additive Gaussian white noise [9-11]. 

The effects of noise on limit cycles can range from frequency shifts [12] to an increased spectral 
bandwidth and a rapid decay of the cycle’s autocorrelation [3, 13]. Quantifying these properties of 
noisy limit cycles is unavoidable if we want to compare the predictions of our models to recorded 
time series. Such a comparison is complicated by the fact that noise can also induce decoherent 
oscillations, in cases where the deterministic model predicts a damped oscillation towards a stable 
fixed point [14-18]. Baxendale and Greenwood [19] show that under certain approximations, such 
quasi-cycles are described by circular orbits whose radius is modulated by an Ornstein-Uhlenbeck 
process [20, 21]. Thompson et al. [22] extend these results to pairs of quasi-cycles sustained by a 
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common noise source. Boland et al. [23] investigate limit cycles of the two-dimensional Brusselator 
system subject to noise, by describing the dynamics along the directions tangential and normal to 
the limit cycle. On a more abstract setting, Teramae et al. [24] and Goldobin et al. [25] derive phase 
equations for noisy limit cycles. 

In the present article, we consider arbitrary autonomous dynamical systems described by an 
ordinary differential equation exhibiting a stable limit cycle, perturbed by additive Gaussian white 
noise. We formally derive an approximation for the dynamics of sample paths that stay close to 
the limit cycle. We formulate these dynamics as a set of coupled stochastic differential equations 
(SDE) for the phase (or longitudinal) coordinate along the limit cycle and a set of complementary 
coordinates collectively referred to as deviation. A similar decomposition was used by Kurrer and 
Schulten [26] to investigate the effects of noise on the Bonhoeffer-van der Pol nonlinear oscillator. 
In analogy to previous work [24], we find that the phase advances at the deterministic speed super¬ 
imposed by a Brownian-motion-like drift. The deviation takes the form of an (n — l)-dimensional 
Ornstein-Uhlenbeck process [20, 21], extending the work of Baxendale and Greenwood [19] on quasi- 
cycles. We apply these results to systems in the normal form of the supercritical Hopf bifurcation. 
The latter describes the emergence of a limit cycle around an unstable focus and appears in numer¬ 
ous ecological and epidemiological models [27-30]. We give explicit approximation formulas for the 
long-term autocovariance and power spectral density of the components of such a process. These 
formulas reproduce the widely observed frequency shift and decoherence of noisy limit cycles, and 
allow a direct comparison of models exhibiting cyclic dynamics to real time series. We asses the 
fidelity of the derived approximations using numerical simulations and exemplify their use for the 
El Nino Southern Oscillation, a well known but poorly understood cyclic climatic phenomenon [31]. 

2 Linear approximation of noisy limit cycles 

Our starting point is a smooth n-dimensional dynamical system ^ = f(y), exhibiting a stable limit 
cycle £. With additive Gaussian white noise, the dynamics take the form of an Ito SDE [32] 

dy = f(y) dt + S dW, (1) 

where W shall be an n-dimensional Wiener process (or standard Brownian motion) with uncorre¬ 
lated components and S G M nxn is some matrix. 

Let any sample path y (t) of the SDE (1) close to £ be split into a longitudinal and a lateral 
component, y (t) = y\ on (t) + yi at (i), where yi on (t) G £ is a point on the limit cycle and yi at (i) is 
perpendicular to the local tangent. Let L(£) be a parameterization of the deterministic limit cycle, 
i.e. such that L(0) = yi O n(0) and = f(L). Let T(f) be the normalized tangent at L(t) in the 
direction of motion and let P{t) be the hyperplane perpendicular to T(f). Abbreviate To = T(0) 
and Po = P(0). The pair (T(i),P(t)) defines a comoving frame along the limit cycle, similar to 
the Frenet frame in 3 dimensions [33]. Let U(t) G SO(n) be an orthogonal matrix, depending 
smoothly on time t, such that T(t) = U(t)To, P(t) = HJ(f)Po and ^U(t)z T P(t) for all z G Po. 
We refer to appendix A for a constructive proof of existence. HJ(0) can be chosen to be the identity 
matrix Id, in which case HJ(t) would be unique. For example, in two dimensions HJ(t) would be the 
rotation that maps To to T(t). See figure 1 for an illustration. Denote V(f) = ^HJ(f). Let P(t) 
and T(t) be the orthogonal projections onto P(t) and the linear span of T(f), respectively. Note 
that T(f), T(t), P(t), P(t), U(t) and V(f) are solely determined by the kinematics of the deterministic 
trajectory L. Let JJ(f) = (Vf)|L( t ) be the Jacobian of the deterministic dynamics around L(t). 

We make the ansatz yi on (f) = L(-r(t)) and yi a t(t) = U(r(f))z(t) for suitable r(f) G M and 
z(t) G Pq. Upon choice of a basis in the (n — l)-dimensional hyperplane Po, z(f) can be described 
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P(t) = U(t )P 0 



Figure 1: Representation of the state y by a phase r and a deviation z £ Pq, illustrated 
for a 3 dimensional phase space. The phase t advances such that the plane perpendicular 
to the tangent T(r), P(r), includes y. The deviation z is defined in the fixed plane Pq. 


by (n — 1) independent variables. Hence, r(t) and z (t) are to be seen as a set of new coordinates 
for y(t) along and perpendicular to the limit cycle, which we shall refer to as phase and deviation, 
respectively. In these coordinates the kinetics formally take the form 


dy dr , .. dr„ r/ . ,dz 

- = - f(L ( T )) + -vHz + u ( r)-, 


( 2 ) 


while we omit to show the explicit dependence on t for brevity. On the other hand, linearising the 
dynamics in the proximity of the limit cycle yields the formal approximation 


^ = f(yion) + J('r)yiat + § + 0{H 11 ylat II 2 ), 

where H is the bound of the second derivative of f (when considered as a bilinear operator I 
M n ), maximized along the entire deterministic limit cycle. Combining (2) with (3) gives 

^ [f(L(r)) + V(r)z] 

=J(t)U(t)z + § —-U(t) —-V(t)z + 0(H 11z11 2 ). 

(JjL lll 


(3) 


(4) 


Equation (4) can be split into the lateral part 




jr(r)U(r)z + § 


dW 

dt 


+ 0(H ||z| 


(5) 


and the longitudinal part 


d T 

--l)[f(L(r))+V(r)z] 


=T(r) 


JJ(r)U(r)z + S 


dW 

dt 


( 6 ) 


— V(t)z + 0(H ||z|| 2 ), 


while we used the fact that V(i)z _L P(t). Equation (6) can be written as 


d-T 

~dt~ l + 


+ o 


(T)U(r)z-¥(r)z + S^,T(r) 
<f(L(r))+V(r)z,T(r)> 


H 11 z I 


(7) 


ll f ( L (T))ll - ||V(r)z| 
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Note that 


||V(t)z 


l|f(L(r))|| 



where s(t ) is the arc-length parameterization along the limit cycle. Furthermore, ||^(r)|| < 
||^(t)|| (see B for a proof). Hence, for small z, (7) can be written as 

dr (J(t)U(t)z-V(t)z + S^,T(t)) 
dt (f(L(r))+¥(r)z,T(r)> 

„ ( H |rf \ <8) 

+ \l|f(L(T))|| 1 - <t(r) z J ' 


where k = is the curvature of the deterministic limit cycle. Summarizing, we can decompose 

sample paths of (1) as 


y = L(r) +U(t)z, 


(9) 


where z and r are correlated stochastic processes whose sample paths approximately solve (5) and 
(8), respectively. 

Due to the local stability of the limit cycle, U T (t)P(t)J(f)U(t) is a stable linear operator on 
Pq. Therefore, (5) describes an Ornstein-Uhlenbeck process on the hyperplane Po, with a possibly 
stochastic noise tensor and Jacobian [34], Intuitively, deviations of sample paths from the limit 
cycle are the result of fluctuations acting against local stabilising dynamics [35], characterized 
predominantly by the limit cycle’s Lyapunov exponent. This is in accordance with results by DeVille 
et al. [36], who showed that for weak noise limit cycles in Hopf normal form have a negative Lyapunov 
exponent. On the other hand, to leading order in z, the phase r advances at the deterministic rate 
modulated by additive white noise, therefore exhibiting a Brownian-motion-like drift away from its 
deterministic value. As such, (8) differs fundamentally from (5), since noise-induced changes of the 
phase are not reversed by the deterministic dynamics. 

Without loss of generality let Po = {(^i, • x n ) £ 1“ : = 0}, so that z = (zo,0) for some 

zo G M n , and let Ho : M n —> M n_1 denote the projection to the first (n — 1) components. Suppose 
that the noise S dW is isotropic with uncorrelated components, i.e. § T § = cr 2 ■ Id for some scalar 
a G M. Then (5) and (8) can be written to leading order as 


dzo ~ Jo(^) z o dt + a dWj, dr 


dt -)- 


a dW p 
l|f(L(T))H’ 


where Jo is defined by 


J 0 (r)z 0 = noU T (r)P(r)J(r)U(r) 



and Wd and W p are Wiener processes of dimension n — 1 and 1, respectively. In fact, Wo and W p 
are uncorrelated, because the two orthogonal projections T(r) and P(r) split the noise S d\V into 
two uncorrelated processes. 


3 Limit cycles emerging through Hopf bifurcations 

In this section, we exemplify the linear approximation (9) for the case where the deterministic limit 
cycle emerges through a supercritical Hopf bifurcation [37, 38]. The ubiquity of studied dynamical 
systems exhibiting a Hopf bifurcation makes this an ideal illustrative example. 
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3.1 Approximation of sample paths 

For simplicity, we consider the normal form [39, §3.4] 



(x 2 + y 2 ) / -Ax/2 - (a - a 0 )y \ (10) 

r 2 \ -Ay/2 + (a - ao)x ) 

+ a dW, 


where a > 0, r > 0 and —A < 0 are the limit cycle’s angular frequency, radius and Lyapunov 
exponent, respectively, ao > 0 is the system’s angular frequency in the proximity of the focus and 
d/0. The first row in (10) describes the linear dynamics in the proximity of the focus, which has 
Lyapunov exponent A/2. The second row describes the nonlinearities giving rise to the limit cycle. 
The Wiener process W appearing in the third row is assumed to have uncorrelated components. In 
polar coordinates ( p , p) the deterministic part of (10) reads 


dp A A p 3 
dt ~ 2 P ~ 2 Z 2 ’ 


dp 

dt 


Ot o + (a 



( 11 ) 


The SDE (10) describes a stationary process with zero mean. Modulo arbitrary phase shifts, the limit 
cycle solution to the deterministic part of (10) is given by L(t) = (r cos at, rsinaf) T . Calculating 
T(f),IP(t), J(t),U(t) and V(f) is straightforward and one eventually obtains from (5), (8) and (9) 
the approximation 


( \ ~ ( ( r + z (t)) ' cos ctr(t) \ 

V y( f ) ) ~ V (r + z(t))- sin ar(t) ) 


for the sample paths, where 


dz 


—Xz dt + <T dWd, 


dr 


m + 2z( ° ~ °°> 

a(r + z) 


dt + 


a dW p 
a(r + z) 


( 12 ) 


(13) 


The SDE for z describes a classical Ornstein-Uhlenbeck process, suggesting that the stationary 
expected squared distance from the deterministic limit cycle will be approximately a 2 /(2A). In the 
following, the ratio NSR = -y/cr 2 /(2A )/r shall be referred to as noise-to-signal ratio (see section 3.2 
for further justification of this name). The 2nd term in the SDE for r vanishes if the focal frequency 
is similar to the limit cycle frequency (ao ~ a), which is a reasonable approximation if the system 
is close to the Hopf bifurcation point. However, for large limit cycle radii this term is expected to 
become non-negligible (see numerical tests in section 4). 


3.2 Leading order approximation of the auto covariance 

We use (12) to derive an approximation for the stationary autocovariance, 

ACV[x]( u) = ^lim E {x(t)x(t + u)} , ( 14 ) 

of the component x. We focus on single components because we wish to draw an analogy to 
ecological times series, which are often only available for a few system variables. We consider the 
SDE (13) for r up to leading order in z. More precisely, we approximate 

x(t) « (r + z) ■ cos (p, (15) 
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where p(t) = at + (a/r)W v is a Brownian motion with deterministic drift and dz = —A z + adW d 
describes a one-dimensional Ornstein-Uhlenbeck process, independent of y>. The process cosy; is 
sometimes referred to as randomised harmonic process [40, §1.2.1]. Inserting (15) into (14) leads to 

r 2 

ACV[x] (u) = lim — E{ cos ip(t) cos ip(t + it)} 

t—t OO 2 k J 

+ lim -E {z(t)z(t + s)} E{ cos <p(t) cos ip(t + it)} 
t—> OO 2 L J 

= ^ [r 2 + ACV[,z](it)] • ACV[cosy>](u), 


where ACV[2;](i4) = is the stationary autocovariance of z [21, §3.3.C]. The autocovariance 

ACV[cosy?](ri) is well known [40, eq. (1.9)], yielding 


ACV[x](u) = — [l + NSR 2 • cos(cm) • e -MW0 2 / 2 . 


(16) 


Here, NSR = yj (<j/r) 2 /(2A) is the noise-to-signal ratio introduced in section 3.1. NSR 2 thus relates 
the signal variance originating in lateral fluctuations around the limit cycle to the variance caused 
by the limit cycle itself. 1/NSR 2 is comparable to the signal-to-noise ratio known from signal 
processing theory [41]. The right-most exponentially decaying factor corresponds to a temporal 
decoherence of the limit cycle [42], due to noise-induced phase drift. This deviation from true 
periodicity in the presence of phase noise, is known as jitter in electronic signal theory [43]. 


3.3 Leading order approximation of the power spectral density 


Similarly to section 3.2, we use the leading order approximation (15) to predict the power spectral 
density, 


PS[x](w) 


of x. Inserting (15) into (17) yields 


1 f T 

lim E — / x(t)e~ lut dt 
T^OO yfT J o 


(17) 


PS[x](w) = A 2 PS[cos<y](w) + PS[^cosy?](o;). 


(18) 


The power spectral density of cosy; is well known [44, §5] and given by 

PS[cosy>](cu) f cos (au)e~^ a ^ r ^ u ^ 2 du 

2 J R 

2 (a/r) 2 [4(a 2 + co 2 ) + (<r/r) 4 ] 

[4(a — co) 2 + (cr/r) 4 ] [4 (a + uj) 2 + (cr/r) 4 ] 

It is straightforward to see that similarly, 

PS[z cos y>](u;) = \ [ ACV[z](u) • cos(au)e~^ a ^ r " >2 ^ u ^ 2 du. 
2 Jr 

Evaluating (20) yields, together with (18) and (19), the power spectral density 

2 r 2 (a/r) 2 [4 (a 2 + uj 2 ) + ( a/r ) 4 ] 

[4(a — co) 2 + (cr/r) 4 ] [4 (a + lo) 2 + (cr/r) 4 ] 


2r 2 ((a/r) 2 + 2A) 

S 

to 

+ uj 2 ) + ((a/r) 2 + 2A) 2 


4(a — co) 2 + ((a/r) 2 + 2A) 2 

4(a + uj) 2 + ((a / r) 2 + 2A) 2 


(19) 


( 20 ) 


( 21 ) 


6 
















The power spectrum (21) can also be obtained from the autocovariance (16) using the Wiener- 
Khintchine theorem [45]. As suggested by (21), the presence of noise leads to a shift of the spectral 
peak to a higher frequency than the limit cycle’s deterministic frequency. Moreover, even for 
very stable limit cycles (A 2> (cr/r) 2 , or NSR ~ 0), the power spectrum retains a non-vanishing 
bandwidth, leading to the temporal decoherence expressed by the decaying autocovariance (16) [46, 
§2.3.1], 

4 Numerical validation 

To test the fidelity of our results, we performed numerical simulations of the exact system (10) and 
its linear approximation (12) over a wide parameter range. We refer to C for technical details. The 
linear approximation is found to have sample paths and a probability distribution that are similar 
to the exact system, provided that noise is sufficiently weak (NSR < 0.1) (figures 2(a,b,c)). For 
stronger noise (NSR > 0.5), the linear approximation has a distribution with heavier tails than 
expected, as well as a stronger peak in the region enclosed by the limit cycle (figures 2(d,e,f)). The 
presence of strong outliers in the linear approximation is due to the linear term — Az in (13), which 
underestimates the force of attraction by the limit cycle in the outer region, compared to the cubic 
term oc — p 3 in the original system (11). The false peak in the inner region is due to the fact that in 
the linear approximation, trajectories with inward deviations exceeding the cycle radius (z < —r ) 
have to pass through the origin (z = —r) on their way back to the limit cycle. This is not the case 
for the exact system, which is attracted to the nearest side of the limit cycle. 

We compared the autocovariances (ACV) and power spectral densities (PS) for the first compo¬ 
nent x of both processes, estimated from generated sample paths. We also compared the computed 
ACV and PS with the leading order formulas (16) and (21), respectively. The linear approximation 
(12), as well as the leading order formulas, reproduce the exact ACV and PS to a great extent, when¬ 
ever (i) the noise is sufficiently weak (NSR < 0.5) and (ii) the focal and limit cycle frequencies are 
similar (|ao — a| < a/10) (figures 3(a,b,d,e)). The approximations fail when a differs significantly 
from ao (|a — ao| > a) and NSR > 0.1 (figures 3(c,f)). This is due to the nonlinear modulation 
of the phase speed at different deviations (rightmost term in (11)), not accurately represented in 
the linear approximation (13). In fact, the leading order formulas were derived by ignoring any 
deviation-dependent modulation of phase speed. Hence, one should expect a reduced accuracy of 
these approximations for systems with strong nonlinearities that modulate their phase speed far 
from the limit cycle. 

5 Conclusions 

The approximation (9), derived for the sample paths of noisy limit cycles, opens the door for a 
qualitative understanding of perturbed systems exhibiting intrinsic cyclic dynamics. These results 
suggest that such systems are better understood by separately considering the dynamics along and 
perpendicular to the deterministic limit cycle, as these are qualitatively different. While deviations 
from the limit cycle decay at a rate determined by the system’s Jacobian, the phase drifts away 
from its deterministic value in a Brownian-motion-like manner. The rate of this random drift is, to 
leading order in the deviation, independent of the cycle’s stability properties. This means that even 
systems with very stable limit cycles, can exhibit low temporal coherence. This becomes clear in 
the approximative formula (16) for the autocovariance of the Hopf normal form, which decays at an 
exponential rate that only depends on the cycle’s radius and the noise strength. This decay makes 
noisy limit cycles fundamentally different from cyclostationary processes, which have a perfect time 
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Figure 2: Sample trajectories computed for the Hopf normal form (10) (left column) and 
its linear approximation (12) (middle column), for weak noise (NSR = 0.1, top row) and 
strong noise (NSR = 0.5, bottom row). The right column shows the probability distribution 
of the x component estimated for both processes. Note that in (c) both curves overlap to 
a great extent. In all cases cto = A = a. The kurtosis [47] in (f) is fa ~ 2.1 for the exact 
process and fa ss 2.6 for its approximation. 


reference and whose autocovariance retains a non-decaying amplitude [51, 52], Therefore, noisy limit 
cycles are a potentially powerful alternative model for temporally decoherent cyclic processes such 
as animal population cycles [53, 54] or the El Nino Southern Oscillation [13]. In fact, the formulas 
for the autocovariance and power spectrum allow, in principle, for a distinction of noisy limit cycles 
from quasi-cycles. The latter are often used to explain observed population cycles [16, 55-57] and 
have well understood autocorrelations and power spectra [16, 58] 

The approximations derived in section 3, both for the sample paths as well as the autocovariance 
and power spectrum, are strictly speaking only valid for systems with isotropic noise in the Hopf 
normal form (10). However, more general systems exhibiting stable limit cycles emerging through 
a Hopf bifurcation are expected to have autocovariances and power spectra that are qualitatively 
similar to the predictions given here, at least in the proximity of the bifurcation point. Hence, 
fitting these formulas to a given cyclic time series might provide a first estimate of the stability of 
the cycle and the amount of perturbations it is subject to. In contrast to conventional nonlinear time 
series autoregressive models [59], the template formulas derived above permit a direct mechanistic 
interpretation of the estimated parameters. This approach is exemplified for a time series of the El 
Nino Southern Oscillation (figure 4(a)), which we hypothesise to be a noisy limit cycle. From the 
fitted autocovariance and power spectrum formulas (figures 4(b,c)), one obtains estimates for the 
Lyapunov exponent and noise variance which can be compared to values obtained from alternative 












Figure 3: Autocovariances (ACV, plots (a,b,c)) and power spectral densities (PS, plots 
(d,e,f)), computed for the Hopf normal form (10) (continuous line) and its linear approxi¬ 
mation (13) (dashed line) using numerical simulations. The formulas (16) for the ACV and 
(21) for the PS are plotted for comparison (dotted line). Note that in figures (a,b,d,e) all 
three curves overlap to a great extent. Parameter values are ao = a and NSR = 0.1 for 
(a,d), «o = a and NSR = 0.5 for (c,e), ao = a/2 and NSR = 0.1 for (c,f). In all cases 
A = a. 


existing models. 

In conclusion, this work provides a starting point for the qualitative understanding and statistical 
validation of stochastic differential equation models, with well-understood deterministic limit cycles. 
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A Lemma on the existence of certain rotations 

Let T (t) E W 1 be a normalized vector that depends smoothly on time t > 0. Let P(t) be the 
hyperplane perpendicular to T(f) and let Pq = P( 0), To = T(0). Then there exists a unique family 
of orthogonal transformations U(£) E SO(n) (t > 0), depending smoothly on time t. such that 


1. U(0) is the identity, 

2. U(f)P 0 = P(t), 

3. U(f)T 0 = T(£), 
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Figure 4: (a) Monthly anomaly of the El Nino Southern Oscillation index N3.4 between 
January 1871 and December 2007 [48]. (b) Estimated autocovariance and (c) periodogram 
estimate of the power spectrum. The dashed curves in (b) and (c) show the formulas for the 
autocovariance (16) and power spectrum (21), respectively, fitted to the data using least 
squares [49]. The ratio er 2 / ACV(0) is estimated at 0.83/yr from the fitted autocovariance 
in (b) and at 0.96/yr from the fitted power spectrum in (c), close to estimates obtained 
from noise-sustained oscillation models (approximately 1.41/yr) by Burgers [13]. The fitted 
cycle period is about 4.2 yr in both cases. The focal Lyapunov exponent A/2 is estimated 
at 0.15/yr from (b) and at 0.17/yr from (c), contrasting previous estimates using delay 
oscillator models (1-1.5/yr) [50]. 


4. ^-(t) z J_ P(t ) for all z € Pq- 


Moreover, U(i) satisfies the linear inhomogeneous differential equation 


dU(t) 

dt 


-T(t)^^U(t)P 0 + 


dm^ T 

~dT L ° 


where Pq is the orthogonal projection onto Pq- 


( 22 ) 


Proof For notational simplicity we will denote X = 4A for any time-dependent variable X. We 
start by showing the existence of HJ(t). Without loss of generality one can assume that 

P 0 = {x = (xi,..,x n ) € M n : xi = 0} = span{e 2 ,..,e n } 


and To = ei, where ei,..,e n is the standard basis in R n . Choose any A (t) 6 SO(n) depending 
smoothly on time and such that A(t)To = T(f) (such a transformation clearly exists). Denote 


S (t) = A T (t)A(t)ei and let 


H (t) 


fSi(t) 

S2(t) 


~S2(t) 

0 


\Sn(t) 0 


where S = (Si ,.., S n ) T . Note that S\(t) = 0, since 


efS(i) = e] C ’A T (t)A(t)ei = —ej r A T (t)A(f)ei 


-S n (t)\ 

0 

0 / 


ej r A T (t)A(t)ei 


-ei’S(t). 


(23) 


In the 2nd step of (23) we used the fact that A T (f)A(i) = — A T (t)A(t), since A (t) is orthogonal. The 
matrix H(t) satishes H T (t) = —H(t), H(f)ei = S(t) and BI(f)Po -L Po- Set D (t) = A(t)W(t)A T (t), 
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then D (t) satisfies D T (i) = —D(i), D(f)T(f) = T(t) and D(f)P(f) _L P(t). Set U(t) as the solution 
to U (t) = D(t)U(t), U(0) = Id, that is, 

HJ(t) = exp f B(s) d-s . 

Uo 

Then T (t) = exp D(.s) ds To = U(f)To, so that claim (3) is satisfied. Moreover, 

U T (f) = exp f D t (s) ds = exp — f D(s) ds = U _1 (f), 

Uo J L 40 

so that HJ(t) is indeed orthogonal. Therefore claim (2) follows from claim (3). Finally, for any z E Pq 
one has U(t)z = D(f)U(f)z _L P(f), by property of D (t). 

We shall now show (22), an immediate consequence of which will be the uniqueness of HJ (t). 
Denote ¥(f) = HJ(t). Due to properties (3) and (4), ¥(f) can be written as 

V(t) = T(f)a T (f) + T(t)Tj (24) 

for some a(f) E Po. Since HJ(f) is orthogonal, we have HJ T (t)¥(i) = —¥ T (t)HJ(f). Hence, for any 
z E Po, 

(a T (f)z)T 0 = (a T (f)z)U T (f)T(f) = U T (f)¥(f)z = -¥ T (f)U(f)z = -T 0 T T (t)U(t)z, 
and thus a = — T T (f)U(t)Po- Therefore (24) can be written as 

¥(t) = -T(t)T T (t)U(t)P 0 + T(t)Tj, 

which proves (22). □ 

B Lemma on the norm of certain rotations 

Let T(f), P(f), To and Po be as in A. Let U (t) E SO(n) depend smoothly on time, such that 
U(f)To = T(f) and ¥(t)z _L P(t) for all z E Po (where ¥ = ^). Then ||¥(t)|| = ||^T(f)||. 

• , 

Proof For notational simplicity we will denote X = -^X for any time-dependent variable X. 
Since To _L Po and T(f) = ¥(i)To, we need to show that ||¥(i)z|| < ||z|| • ||¥(f)To|| for all z E Po- 
Since ¥(f)z and T(f) are parallel, one has 

l|¥(t)z|| = |(¥(f)z, T(f))| = |(¥(f)z,U(f)T 0 )| = |(U T (f)¥(f)z, T 0 )| . 

Since HJ (t) is orthogonal, one has U T (i)¥(t) = ¥ T (f)U(f). Hence 

||¥(t)z|| = |(U(f)z,¥(f)T 0 )| < ||U(f)z|| • ||¥(f)T 0 || = ||z|| • ||¥(f)T 0 ||, 

as claimed. □ 
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C Details on numerical methods 


For the generation of sample paths, we used an explicit two-step Runge-Kutta scheme of mean 
square order 3/2 [60, §3.4, Theorem 3.3], implemented in C++. We normalized r = 1 and a = 2 tt, 
and considered ao £ [0.1,10] • a and A G [0.1,10] • a. Noise to signal ratios were considered in the 
range NSR G [0.001,1]. The integration time step was set to 10 -4 cycle periods; decreasing it did 
not significantly change the outcome of the simulations. Autocovariances were estimated from the 
sample autocovariances of generated paths, spanning 10 6 points over 10 4 cycle periods [61, §2.5.2], 
Power spectral densities were estimated by averaging the sample periodograms of 100 independent 
sample paths, each spanning 10 4 points over 10 2 cycle periods [61, §13.1]. Probability distributions 
were estimated from 10 5 points spanning 10 3 cycle periods using a Gaussian kernel density estimator 
and Silverman’s rule of thumb [62, p. 48, eq. (3.31)]. 
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